Dissipative particle dynamics with energy 

conservation 



J.Bonet Avalos and A.D. Mackie 
Departament d'Enginyeria Qmmica, ETSEQ, Universitat Rovira i Virgili. 
Carretera de Salou s/n, 43006 Tarragona (Spain) 

February 1, 2008 

Abstract 

The stochastic differential equations for a model of dissipative particle dynamics 
with both total energy and total momentum conservation in the particle-particle 
interactions are presented. The corresponding Fokker-Planck equation for the evo- 
lution of the probability distribution for the system is deduced together with the 
corresponding fluctuation-dissipation theorems ensuring that the ab initio chosen 
equilibrium probability distribution for the relevant variables is a stationary solu- 
tion. When energy conservation is included, the system can sustain temperature 
gradients and heat flow can be modeled. 

PACS 47.11+j -Computational methods in fluid dynamics 

PACS 05.70Ln -Nonequilibrium thermodynamics, irreversible processes 

PACS 02.70Ns -Molecular dynamics and particle methods 

Much attention has recently been paid to the simulation strategy for the dynamics 
of complex fluids known as Dissipative Particle Dynamics (DPD), which was first intro- 
duced by Hoogerbrugge et In this model, a system of mesoscopic particles can 
interact via direct conservative potentials, as in molecular dynamics (MD) simulations 
but, in addition, the particles exert friction and brownian forces on each other. The dissi- 
pative and random interactions are chosen in such a way that the center of mass motion 
of each interacting pair is insensitive to their action. In this way, the system relaxes to 
its equilibrium much faster than in MD simulations and, at the same time, the interac- 
tion conserves the total momentum. This second feature allows the system to exhibit a 
hydrodynamic behaviour from a macroscopic point of view. The model is isotropic and 
Galilean- invariant, in contrast with models denned on a lattice, and has the potential to 
be computationally efficient. 

Thus, DPD appears as an interesting tool although its capability for making quanti- 
tative predictions on the dynamics of complex systems remains still to be explored. Some 
steps forward, however, have already been taken. For instance, Espanol et have 



derived the fTuctuation-disipation theorem appropriate to these dissipative particles, the 
true hydrodynamic behaviour of the model has been proved 0] and the relationship of 
the macroscopic transport coefficients with the mesoscopic parameters defining the model 
has been established 0]. The DPD method has also been applied with remarkable suc- 
cess to account for the many-body problem of the hydrodynamic interactions in colloidal 
suspensions as well as to domain growth in binary immiscible fluids Q. 

As originally formulated, DPD can only deal with isothermal conditions since the 
energy is not conserved in the interaction 0]. The system therefore cannot sustain tem- 
perature gradients and, hence, no heat flow can be modeled^]. This limitation excludes 
the DPD simulation strategy from those problems in which non-equilibrium situations 
involving temperature gradients or heat generation are important, as is the case for con- 
vection or in chemical reactions. The aim of this letter is to extend the DPD algorithm to 
account for both momentum and energy conservation in the particle-particle interaction, 
whilst maintaining the irreversible nature of these interactions. 

The dissipative particles can be viewed as /umps|Q or clusters of true physical particles 
or as particles with internal structure bearing some degrees of freedom. The DPD model 
is mesoscopic in nature since it resolves only the overall center-of-mass motion of the 
cluster and avoids the description of the variables specifying its internal state. We will 
account for the energy stored in the internal degrees of freedom of the particles without 
explicit consideration of any internal Hamiltonian, in a model inspired by the treatment 
of hydrodynamic fluctuations |7|, |], ||. 

Our model is based on the following assumptions: 



1. The system contains iV particles interacting with each other via conservative as 
well as dissipative interactions. The conservative interactions are described by the 
Hamiltonian 

N ( 2 N ") 

H({ri}, {Pi}) = E | + 5>(ry) + r*\n) | (1) 

where = |r*j — fj\. The Hamiltonian depends on the momenta and positions of all 
the particles. The particles interact with each other through pair potentials, ip{ r ij), 
which depend only on the distance between the particles, and with an external field, 

2. In addition, the particles can store energy in some internal degrees of freedom. The 
internal energy, Ui with m > 0, is introduced as a new relevant coordinate. Together 
with the momentum pi and the position r*j this completely specifies the state of the 
dissipative particle at a given instant t. 

3. The particle-particle interaction is pairwise and conserves the total momentum, the 
total angular momentum and the total energy when the internal energy of the pair 
is taken into account. 

4. The internal states of the particle have no dynamics in the sense that they are 
always in equilibrium with themselves. This allows us to define a function Sj(wj). 



2 



This function can arbitrarily be chosen according to the user's need, only constrained 
to the requirements: a) must be a differentiable monotonously increasing function 
of so that the function Ui(si) exists and 9i = dui/dsi exists and is always positive, 
and b) Si(uj) is a concave function of its argument |10||. Defined in this way, Sj can 



be viewed as a mesoscopic entropy of the i th particle, and 9i can be seen as the 
particle's temperature. The change in Ui and in Sj are related by a Gibbs equation]?] 
Oidsi = dui, which implies 9is\ = Ui, where the dot over the variables is used to 
denote time-differentiation from now on. 

5. In the absence of random terms, the pairwise particle-particle interaction is irre- 
versible and must satisfy Sj + Sj > 0, where % and j label an arbitrary pair of 
particles in interaction with each other. 

6. Ui, Si and 9% must remain unchanged under a Galilean transformation, so that these 
variables are true scalars. 

7. The equilibrium probability distribution for the relevant variables of the system 
under isothermal conditions is chosen to be 

N 

Pedrt}, {Pi}, {ui}) ~ e -Hmhm)/Kr rj ^ma-^at ^ 

i=i 

where k is Boltzmann's constant and T is the thermodynamic temperature. The first 
factor on the right hand side of eq. (0) contains the probability distribution for the 
variables {pi} as given from equilibrium statistical mechanics. The second fac- 
tor on the right hand side corresponds, in turn, to the probability distribution for the 
internal energy of the particles as obtained from equilibrium fluctuation theory fT0| . 



The maximum of Si(ui)/k — Ui/kT takes place at 9i = T, in agreement with our 
interpretation of 9i as the particle's temperature. Once the equilibrium probability 
distribution is set, all the thermodynamic properties of the model can be derived 
from the partition function Z(T, V, N) = f(J\i dpidridui) P e ({rj}, {p>i}, 



We proceed to the analysis of the dynamics of the model defined so far, although 
the details of the calculation will be given elsewhere. For simplicity, we will analyse 
two particles, i and j say, due to the pairwise additivity of the interaction, and give the 
complete expression at the end. Since the state of the system is specified by the set 
{pi} and {ui}, we have to supply equations for the evolution of these variables. The 
change in the position and momentum of the i th particle due to the interactions with the 
j th particle is given by 

rl = ^ (3) 
m 

'p % = fV + FT' + Fg + Fg (4) 

where F-j = —dip^r^/dfi and Ff xt = —dip ext {fi)/dri are the forces due to the conserva- 
tive interactions. F£j is, by construction, directed along the unit vector f,y = (r} — r^/r^. 
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Fg stands for the dissipative particle-particle interaction force and Fg is the random 
force associated with the former. The mechanisms driving the change in the internal 
energy of the particles are assumed to be of two kinds. On the one hand, the work done 
by the dissipative forces increases the internal energy of the interacting particles. Since 
they are identical, we assume that this work is shared in equal amounts by the particles 
and is irrespective of the temperatures 0$ and 9j. At the same time, the action of the 
random force cools the particles transfering internal energy back to mechanical energy. 
We want, in addition, that viscosity and thermal conductivity can be independently mod- 
eled. Hence, we also consider, on the other hand, that interacting particles can exchange 
internal energy if 9i 7^ 9j, something that we can call mesoscopic heat flow qg between 
particles. Associated with this dissipative current, we also add a random heat flow qg. 
Thus, we write 

Ui = ^ipj - Pi) ■ (Fg + Fg) + qg + qg (5) 

The dissipative and random terms have to be of the form Fg' R ~ r y - and qg' R = —qfi' R - 
Otherwise, the requirements of point 3 would be violated. 

The analysis of point 5, permits the identification of the so-called thermodynamic 
/orcesQ causing the dissipative currents to exist. In our case, it is rather intuitive that 
the friction forces are due to the momentum difference between interacting particles, and 
that the heat flow is due to a temperature difference. We, thus, assume linear relationships 
of the form 

Fg = C(r lj )^r ij r ij ■ (pj - p { ) and qg = A(r y ) f ^ - ^-J (6) 

where C( r ij) an d A(r y -) are arbitrary positive coefficients, that are even functions under 
time-reversal ||1 1|| . A convenient choice is to assume them to be functions of the slow 



variable only. Galilean invariance of Fg and qg is thus guaranteed. In this way, all 
the deterministic interactions are completely specified. 

The properties of the random terms are chosen to also parallel the theory of hydro- 
dynamic fluctuations m [!|. Since Fg and qg are not coupled, we will demand that the 
random terms Fg and qg be statistically independent. They can be written in the form 

'■g '■;,]-,, :F t jit ) and $ = Sign(i - j) Aij Qijit) (7) 
where the function Sign(i — j) is 1 if i > j and — 1 if i < j, ensuring that g^ = —q R . The 



scalar random variables T%j and Qij are stationary, Gaussian and white [|T2|, |TT|j , with zero 
mean and correlations 

{^{t)W)) = (Q l3 {t)Q kl {t')) = (5 ik 5 jt + 5 u S jk ) 5(t - t') (8) 

Tij and Ay are functions to be determined later. Note that ( and T^- are, respectively, 
7a>£> and aujR in ref.0. 

The stochastic differential equations eqs. (|3|), (^) and (|5|), together with eqs. (Q) and 
the properties of the random terms given in eqs. ([?]) and lead to the Fokker-Planck 
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equation 



§^ p ({nh {£}, {«*}) = £ core P({n}, {«*}) + ^^({n}, {£}, M) (9) 

where the convective operator, L con , and the diffusive operator, L dl ^ , are defined as 
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where we have in addition defined the operator 



d d 



+ 



(10) 




+ ^®-pi 




(12) 



Imposing that the equilibrium distribution function given in eq. (g) is a stationary solution 
of eq. @ Q, we derive the corresponding fluctuation-dissipation theorems 



r? 



Ik &ij Cij 



K 3 = 2* Ay 



(13) 
(14) 



where we have defined the mean inverse temperature as 0" 1 = (1/0, + l/6 , J )/2, that is 
a function only of Ui and Uj. Eq. (|9|) together with eqs. ( JT3] ) and ( |T4] ) are a generaliza- 
tion of the results found in ref.0, which refer to isothermal conditions and momentum 
conservation only. A crucial difference regarding the momentum change is, however, the 
fact that the fluctuation-dissipation theorem given in eq. ( |13|) relates the strength of 
the random force with the temperature 0^-, instead of the thermodynamic temperature 
T0. Thus, the model presented here is defined in terms of particle properties only, with 
no reference to macroscopic magnitudes such as the thermodynamic temperature or the 
density of the system. This neat property permits the use of DPD in situations other 
than isothermal. Together with the updating algorithms shown below and the second 
fluctuation-dissipation, eq. (|I4]), this is the main result of this letter. 



Note that the Langevin equations eqs. 



and 



are subject to the so-called Ito- 
Stratonovich dilemma |L3] due to the occurence of products of fast variables, such as pi or 
0y , and random variables which are 5-correlated in time. We take, however, eq. ([]) as the 
true definition of the random processes. The choice made here for these random processes 



1 In fact, the probability flux must also vanish for a system in thermodynamic equilibrium JlT[| 
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is based on the detailed balance condition and on the weak coupling assumpUon\\ll\\, which 
permit that the properties of the random forces be independent of overall thermodynamic 
properties. Finally, the updating algorithms can directly be obtained from eq. (|9|), using 
eqs. ( |T3|) and (14), with no ambiguity After some algebra, we obtain the new values of 
the variables r ■ = f{(t + 6t), p \ = Pi(t + St) and u\ = Ui(t + 5t) in terms of the old ones 
at t 



P 



Ti H St 

m 

{ 3& 



(15) 



Fy + (Pj ~ Pi) ■ fijfij 
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(17) 



where we have defined an effective friction coefficient 

k ( d d 



i + 



+ 



6; 



2 \dui ' duj J 13 



toff and n^' are independent random numbers defined from the random variables J-^ 



ij 



and Qij and are Gaussian with zero mean and variance (f^Qj^) = (fi^O^j ) = (5^5 

0~il°~jk)- 



The model presented in this letter constitutes a generalization of the DPD method. 
The addition of energy conservation in the particle-particle interaction in a consistent way, 
allows the derivation of an updating algorithm defined in terms of particle variables only. 
We have obtained two fluctuation-dissipation theorems which permit that the proper ther- 
modynamic equilibrium be reached for the model under appropriate conditions. It should 
be mentioned that, with respect to previous treatments 0, the fluctuation-dissipation 
theorem for the random force derived in this work contains the local temperature 
instead of the thermodynamic temperature T, stressing the local nature of our model. 
In addition, the fluctuation-dissipation theorem for the heat flux has no counterpart in 
previous DPD formulations. The corresponding computer implementation of the model 
has already shown that the total energy is well conserved for an isolated system, although 
some dependency on the time step St has been observed. Therefore, our model can sustain 
temperature gradients and thus heat flow can be described. 
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